Enrichment x stress MA
Supplementary Material
Setting-up
Loading packages
# devtools::install_github('Mikata-Project/ggthemr', force = TRUE)
pacman::p_load(tidyverse, here, metafor, clubSandwich, orchaRd, MuMIn, patchwork,
GoodmanKruskal, networkD3, ggplot2, visdat, ggalluvial, ggthemr, cowplot, grDevices,
png, grid, gridGraphics, pander, formatR)
# needed for model selection using MuMIn within metafor
eval(metafor:::.MuMIn)Loading data and functions
This loads the unprocessed datafile and custom functions including
- calculating ‘focal’ and ‘pair-wise’ effect sizes and variance
- calculating SMD
- creating VCV of variance
dat <- read_csv(here("Data", "Data_raw.csv"))
# Load custom function to extract data
source(here("R/Functions.R"))Data organisation
- removing study (Wang et al_2020) with negative values as lnRR cannot be calculated with negative values
- rounding down sample sizes that are reported as decimals due to averaging n across treatments
- getting effect sizes from function
- ‘flipping’ effect sizes so that all effect sizes are higher values = individuals do better in cognitive assays
- assigning human readable terms, and creating VCV of variance
# removing study with negative values as these are unable to be used for lnRR
dat <- droplevels(dat[!dat$First_author == 'Wang',])
#rounding down sample sizes
dat$CC_n <- floor(dat$CC_n)
dat$EC_n <- floor(dat$EC_n)
dat$CS_n <- floor(dat$CS_n)
dat$ES_n <- floor(dat$CS_n)
# 'Focal' effect_size
effect_size <- with(dat, mapply(effect_set,
CC_n ,
CC_mean,
CC_SD,
EC_n,
EC_mean,
EC_SD,
CS_n,
CS_mean,
CS_SD,
ES_n,
ES_mean,
ES_SD,
percent = Response_percent,
SIMPLIFY = FALSE))
effect_size <- map_dfr(effect_size, I)
# 'Pairwise' effect size
effect_size2 <- with(dat, mapply(effect_set2,
CC_n ,
CC_mean,
CC_SD,
EC_n,
EC_mean,
EC_SD,
CS_n,
CS_mean,
CS_SD,
ES_n,
ES_mean,
ES_SD,
percent = Response_percent,
SIMPLIFY = FALSE))
effect_size2 <- map_dfr(effect_size2, I)
full_info <- which(complete.cases(effect_size) == TRUE)
# adding effect sizes as column
dat <- bind_cols(dat, effect_size, effect_size2)
dat <- dat[full_info, ]
#Flipping 'lower is better' to 'higher is better' effect sizes
#flipping lnRR for values where higher = worse
dat$lnRR_Ea <- ifelse(dat$Response_direction == 2, dat$lnRR_E*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E))
# currently NAswhich causes error
dat$lnRR_Sa <- ifelse(dat$Response_direction == 2, dat$lnRR_S*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S)) # currently NAswhich causes error
dat$lnRR_ESa <- ifelse(dat$Response_direction == 2, dat$lnRR_ES*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_ES)) # currently NAswhich causes error
#flipping 'pure effect sizes'
dat$lnRR_E2a <- ifelse(dat$Response_direction == 2, dat$lnRR_E2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E2)) # currently NAswhich causes error
dat$lnRR_S2a <- ifelse(dat$Response_direction == 2, dat$lnRR_S2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S2)) # currently NAswhich causes error
dat$lnRR_ES2a <- ifelse(dat$Response_direction == 2, dat$lnRR_ES2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_ES2)) # currently NAswhich causes error
dat$lnRR_E3a <- ifelse(dat$Response_direction == 2, dat$lnRR_E3*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E3)) # currently NAswhich causes error
dat$lnRR_S3a <- ifelse(dat$Response_direction == 2, dat$lnRR_S3*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S3)) # currently NAswhich causes error
#flipping SMD
dat$SMD_Ea <- ifelse(dat$Response_direction == 2, dat$SMD_E*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_E)) # currently NAswhich causes error
dat$SMD_Sa <- ifelse(dat$Response_direction == 2, dat$SMD_S*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_S)) # currently NAswhich causes error
dat$SMD_ESa <- ifelse(dat$Response_direction == 2, dat$SMD_ES*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_ES))
# assigning human readable terms
dat <- dat %>% mutate(Type_assay = case_when(Type_assay == 1 ~ "Habituation",
Type_assay == 2 ~ "Conditioning",
Type_assay == 3 ~ "Recognition",
Type_assay == 4 ~ "Unclear"),
Learning_vs_memory = case_when(Learning_vs_memory == 1 ~ "Learning",
Learning_vs_memory == 2 ~ "Memory",
Learning_vs_memory == 3 ~ "Habituation"),
Type_reinforcement = case_when(Type_reinforcement== 1 ~"Appetitive",
Type_reinforcement== 2 ~ "Aversive",
Type_reinforcement== 3 ~ "Not applicable",
Type_reinforcement== 4 ~ "Unclear"),
Type_stress_exposure = case_when(Type_stress_exposure == 1 ~ "Density",
Type_stress_exposure == 2 ~ "Scent",
Type_stress_exposure == 3 ~ "Shock",
Type_stress_exposure == 4 ~ "Exertion",
Type_stress_exposure == 5 ~ "Restraint",
Type_stress_exposure == 6 ~ "MS",
Type_stress_exposure == 7 ~ "Circadian rhythm",
Type_stress_exposure == 8 ~ "Noise",
Type_stress_exposure == 9 ~ "Other",
Type_stress_exposure == 10 ~ "Combination",
Type_stress_exposure == 11 ~ "unclear"),
Age_stress_exposure = case_when(Age_stress_exposure == 1 ~ "Prenatal",
Age_stress_exposure == 2 ~ "Early postnatal",
Age_stress_exposure == 3 ~ "Adolescent",
Age_stress_exposure == 4 ~ "Adult",
Age_stress_exposure == 5 ~ "Unclear"),
Stress_duration = case_when(Stress_duration == 1 ~ "Acute",
Stress_duration == 2 ~ "Chronic",
Stress_duration == 3 ~ "Intermittent",
Stress_duration == 4 ~ "Unclear"),
EE_social = case_when(EE_social == 1 ~ "Social",
EE_social== 2 ~ "Non-social",
EE_social == 3 ~ "Unclear"),
EE_exercise = case_when(EE_exercise == 1 ~ "Exercise",
EE_exercise == 2 ~ "No exercise"),
Age_EE_exposure = case_when(Age_EE_exposure == 1 ~ "Prenatal",
Age_EE_exposure == 2 ~ "Early postnatal",
Age_EE_exposure == 3 ~ "Adolescent",
Age_EE_exposure == 4 ~ "Adult",
Age_EE_exposure == 5 ~ "Unclear"),
Exposure_order = case_when(Exposure_order == 1 ~ "Stress first",
Exposure_order == 2 ~ "Enrichment first",
Exposure_order == 3 ~ "Concurrently",
Exposure_order == 4 ~ "Unclear"),
Age_assay = case_when(Age_assay == 1 ~ "Early postnatal",
Age_assay == 2 ~ "Adolescent",
Age_assay == 3 ~ "Adult",
Age_assay == 4 ~ "Unclear"),
Sex = case_when(Sex == 1 ~ "Female",
Sex == 2 ~ "Male",
Sex == 3 ~ "Mixed",
Sex == 4 ~ "Unclear"),
Type_EE_exposure = case_when(Type_EE_exposure == 1 ~ "Nesting material",
Type_EE_exposure == 2 ~ "Objects",
Type_EE_exposure == 3 ~ "Cage complexity",
Type_EE_exposure == 4 ~ "Wheel/trademill",
Type_EE_exposure == 5 ~ "Combination",
Type_EE_exposure == 6 ~ "Other",
Type_EE_exposure == 7 ~ "Unclear"),
ROB_blinding = case_when(ROB_blinding == 1 ~ "Yes",
ROB_blinding == 2 ~ "No",
ROB_blinding == 3 ~ "Unclear"),
ROB_randomisation = case_when(ROB_randomisation == 1 ~ "Yes",
ROB_randomisation == 2 ~ "No",
ROB_randomisation == 3 ~ "Unclear"))
#making variance VCVs
VCV_E <- impute_covariance_matrix(vi = dat$lnRRV_E, cluster = dat$Study_ID, r = 0.5)
VCV_S <- impute_covariance_matrix(vi = dat$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
VCV_ES <- impute_covariance_matrix(vi = dat$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)
VCV_Ea <- impute_covariance_matrix(vi = dat$SMDV_E, cluster = dat$Study_ID, r = 0.5)
VCV_Sa <- impute_covariance_matrix(vi = dat$SMDV_S, cluster = dat$Study_ID, r = 0.5)
VCV_ESa <- impute_covariance_matrix(vi = dat$SMDV_ES, cluster = dat$Study_ID, r = 0.5)
#write.csv(dat, file = here("Data", 'Data_processed.csv'), row.names = TRUE)Data exploration
General
Number of effect sizes
# Number of effect sizes
length(unique(dat$ES_ID))## [1] 92
Number of studies
length(unique(dat$Study_ID))## [1] 30
Publication year range
min(dat$Year_published)## [1] 2006
max(dat$Year_published)## [1] 2021
Visual of missing data
plot_missing <- vis_miss(dat) + theme(plot.title = element_text(hjust = 0.5, vjust = 3),
plot.margin = margin(t = 0.5, r = 3, b = 1, l = 1, unit = "cm")) + ggtitle("Missing data in the selected predictors") #no missing values
plot_missing# useGoodman and Kruskal’s τ measure of association between categorical
# predictor variables (function from package GoodmanKruskal:
# https://cran.r-project.org/web/packages/GoodmanKruskal/vignettes/GoodmanKruskal.html)
# GKmatrix <- GKtauDataframe(subset(dat, select = c('Sex', 'Type_assay',
# 'Learning_vs_memory', #'Type_reinforcement', 'Type_stress_exposure',
# 'Age_stress_exposure', 'Stress_duration', #'EE_social_HR', 'EE_exercise',
# 'Age_EE_exposure', 'Exposure_order', 'Age_assay'))) plot(GKmatrix)
# simple pairwise contingency tables table(dat$Type_assay,
# dat$Type_reinforcement) table(dat$Age_stress_exposure, dat$Age_EE_exposure)
# table(dat$Type_stress_exposure, dat$Age_stress_exposure)
# table(dat$Type_stress_exposure, dat$Age_assay)
# table(dat$Type_stress_exposure, dat$Stress_duration)Visual of missing data for each column. Note that most missing data are for different descriptive statistics as papers are often report results using different types of descriptive statistics (i.e., SD, SE, mean, median etc). Cannot reduce text size in this figure
Alluvial diagrams
Shows the relationship/nestedness of different data elements (used to produce Fig. 2)
Subjects info: species-strain-sex
freq_A <- as.data.frame(table(dat$Sex, dat$Common_species, dat$Strain)) %>%
rename(Sex = Var1, Species = Var2, Strain = Var3) #make a data frame of frequencies for three selected variables
is_alluvia_form(as.data.frame(freq_A), axes = 1:3, silent = TRUE)
p1 <- ggplot(data = freq_A, aes(axis1 = Sex, axis2 = Species, axis3 = Strain, y = Freq)) +
geom_alluvium(aes(fill = Sex)) + geom_flow() + geom_stratum(aes(fill = Sex)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + geom_text(stat
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + "stratum",
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + aes(label
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + after_stat(stratum)))
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + #theme_minimal()
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
theme_void() + theme(legend.position = "none", plot.title = element_text(hjust = 0,
vjust = 3), axis.title.x = element_text(), axis.text.x = element_text(face = "bold"),
plot.margin = unit(c(1, 1, 0, 1), "cm")) + scale_x_discrete(limits = c("Sex",
"Species", "Strain"), expand = c(0.15, 0.05), position = "top") + ggtitle("A study subjects")
p1 Fig. 2a Relationship/nestedness between study species, strain of rodent, and sex
Stress info: age-duration-type stress
freq_C <- as.data.frame(table(dat$Age_stress_exposure, dat$Stress_duration, dat$Type_stress_exposure)) %>%
rename(Age_stress = Var1, Duration_stress = Var2, Type_stress = Var3) #make a data frame of frequencies for three selected variables
is_alluvia_form(as.data.frame(freq_C), axes = 1:3, silent = TRUE)
p3 <- ggplot(data = freq_C, aes(axis1 = Age_stress, axis2 = Duration_stress, axis3 = Type_stress,
y = Freq)) + geom_alluvium(aes(fill = Age_stress)) + geom_flow() + geom_stratum(aes(fill = Age_stress)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + geom_text(stat
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + "stratum",
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + aes(label
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + after_stat(stratum)))
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + #theme_minimal()
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
theme_void() + theme(legend.position = "none", plot.title = element_text(hjust = 0,
vjust = 3), axis.title.x = element_text(), axis.text.x = element_text(face = "bold"),
plot.margin = unit(c(1, 1, 0, 1), "cm")) + scale_x_discrete(limits = c("Age",
"Duration", "Type"), expand = c(0.1, 0.1), position = "top") + ggtitle("C stress exposure")
p3 Fig. 2c Relationship/nestedess between the age of stress exposure, the duration of the stres, and the type of stress
Assay info: LvsM-type-reinforcement
freq_D <- as.data.frame(table(dat$Learning_vs_memory, dat$Type_assay, dat$Type_reinforcement)) %>%
rename(Learning_Memory = Var1, Type = Var2, Reinforcement = Var3) #make a data frame of frequencies for three selected variables
is_alluvia_form(as.data.frame(freq_D), axes = 1:3, silent = TRUE)
p4 <- ggplot(data = freq_D, aes(axis1 = Learning_Memory, axis2 = Type, axis3 = Reinforcement,
y = Freq)) + geom_alluvium(aes(fill = Learning_Memory)) + geom_flow() + geom_stratum(aes(fill = Learning_Memory)) +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + geom_text(stat
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + "stratum",
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + aes(label
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + =
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + after_stat(stratum)))
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + #theme_minimal()
geom_text(stat = "stratum", aes(label = after_stat(stratum))) + #theme_minimal() + +
theme_void() + theme(legend.position = "none", plot.title = element_text(hjust = 0,
vjust = 3), axis.title.x = element_text(), axis.text.x = element_text(face = "bold"),
plot.margin = unit(c(1, 1, 0, 1), "cm")) + scale_x_discrete(limits = c("Learning_Memory",
"Type", "Reinforcement"), expand = c(0.1, 0.1), position = "top") + ggtitle("D cognitive assay")
p4 Fig. 2d Relationship/nestedness between if the assay measured a learning or memory trait, the type of assay used, and the type of reinforcement used
Combined plot
# p1 + scale_fill_brewer(palette = 'Set3') #Pastel1
(p1 + scale_fill_brewer(palette = "Set3"))/(p2 + scale_fill_brewer(palette = "Set3"))/(p3 +
scale_fill_brewer(palette = "Set3"))/(p4 + scale_fill_brewer(palette = "Set3")) +
plot_layout(ncol = 1, heights = c(1, 1, 1, 1, 1))# ggsave(file = './figs/Alluvial_diagrams_v2.pdf', width = 10, height = 12,
# units = 'cm', dpi = 300, scale = 2) #, device = cairo_pdf)Fig. 2 The combined plot of all alluvial diagrams used to make Fig. 2 in the main text
Risk of Bias
Percent of studies that used randomisation
dat %>%
group_by(ROB_randomisation) %>%
summarise(n = n_distinct(Study_ID))## # A tibble: 2 x 2
## ROB_randomisation n
## <chr> <int>
## 1 Unclear 15
## 2 Yes 15
15/30## [1] 0.5
Percent of studies that used blinding
# blinding
dat %>%
group_by(ROB_blinding) %>%
summarise(n = n_distinct(Study_ID))## # A tibble: 2 x 2
## ROB_blinding n
## <chr> <int>
## 1 Unclear 24
## 2 Yes 6
6/30## [1] 0.2
Modelling with lnRR
Main effect: environmental enrichment
Meta-analysis
#dat <- read_csv(here("Data","Data_processed.csv"))
mod_E0 <- rma.mv(yi = lnRR_Ea, V = VCV_E, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.7836 19.5672 27.5672 37.6106 28.0323
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0082 0.0905 30 no Study_ID
## sigma^2.2 0.0345 0.1858 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 809.2712, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1860 0.0320 5.8116 91 <.0001 0.1224 0.2496 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.310673e-01 1.786605e-01 7.524068e-01 6.097041e-10
orchard_plot(mod_E0, mod = "Int", xlab = "lnRR", alpha=0.4) + # Orchard plot
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5)+ # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2)+ # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_colour_manual(values = "darkorange")+ # change colours
scale_fill_manual(values="darkorange")+
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) Fig. 3 Orchard plot showing meta-analytic mean and 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Meta-regression: uni-moderator
Learning vs Memory
Was the response learning or memory?
mod_E1 <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E1) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.3793 18.7586 28.7586 41.2577 29.4729
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0076 0.0873 30 no Study_ID
## sigma^2.2 0.0340 0.1843 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 802.5794, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 18.2861, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning 0.2303 0.0450 5.1227 90 <.0001 0.1410
## Learning_vs_memoryMemory 0.1624 0.0355 4.5713 90 <.0001 0.0919
## ci.ub
## Learning_vs_memoryLearning 0.3197 ***
## Learning_vs_memoryMemory 0.2330 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E1) ## R2_marginal R2_coditional
## 0.02572583 1.00000000
LvsM_E<- orchard_plot(mod_E1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21) + # mean estimate
scale_size_continuous(range = c(1, 7)) + # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
LvsM_E Fig. 4a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparison
contra_mod_E1 <- contrast_fun(data = dat, response = lnRR_Ea, moderator = Learning_vs_memory,
VCV = VCV_E)
res_table_mod_E1 <- get_estimate(model = mod_E1, contra = contra_mod_E1, moderator = Type_assay)
res_table_mod_E1## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Learning_vs_memoryLearn… 0.230 0.141 0.320 1.70e-6 -0.185 0.645
## 2 Learning_vs_memoryMemory 0.162 0.0919 0.233 1.54e-5 -0.249 0.574
## 3 Learning_vs_memoryLearn… -0.0679 -0.164 0.0284 1.65e-1 -0.484 0.349
Type of assay
The broad category of the type of assay used to measure learning or memory
dat1 <- filter(dat, Type_assay %in% c("Recognition", "Habituation", "Conditioning"))
VCV_E1 <- impute_covariance_matrix(vi = dat1$lnRRV_E, cluster = dat$Study_ID, r = 0.5)
mod_E2 <- rma.mv(yi = lnRR_Ea, V = VCV_E1, mod = ~Type_assay-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_E2)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -7.8496 15.6993 27.6993 42.6311 28.7237
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0099 0.0995 30 no Study_ID
## sigma^2.2 0.0321 0.1792 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 661.8903, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.7993, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Type_assayConditioning 0.2167 0.0351 6.1650 89 <.0001 0.1468 0.2865
## Type_assayHabituation 0.1821 0.1147 1.5886 89 0.1157 -0.0457 0.4100
## Type_assayRecognition 0.0554 0.0640 0.8659 89 0.3889 -0.0717 0.1826
##
## Type_assayConditioning ***
## Type_assayHabituation
## Type_assayRecognition
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E2) ## R2_marginal R2_coditional
## 0.07376134 1.00000000
Learning_E <- orchard_plot(mod_E2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34", "grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Learning_E Fig. 4b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparisons
contra_mod_E2 <- contrast_fun(data = dat1, response = lnRR_Ea, moderator = Type_assay,
VCV = VCV_E1)
res_table_mod_E2 <- get_estimate(model = mod_E2, contra = contra_mod_E2, moderator = Type_assay)
res_table_mod_E2## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Conditioning 0.217 0.147 0.287 2.02e-8 -0.197 0.630
## 2 Habituation 0.182 -0.0457 0.410 1.16e-1 -0.285 0.649
## 3 Recognition 0.0554 -0.0717 0.183 3.89e-1 -0.371 0.482
## 4 Conditioning-Habitua… -0.0345 -0.262 0.193 7.63e-1 -0.501 0.432
## 5 Conditioning-Recogni… -0.161 -0.296 -0.0265 1.96e-2 -0.590 0.268
## 6 Habituation-Recognit… -0.127 -0.381 0.127 3.24e-1 -0.607 0.353
Type of reinforcement
If conditioning was used, was aversive or appetitive reinforcement used?
dat2 <- filter(dat, Type_reinforcement %in% c("Appetitive", "Aversive", "Not applicable"))
VCV_E2 <- impute_covariance_matrix(vi = dat2$lnRRV_E, cluster = dat2$Study_ID, r = 0.5)
mod_E3 <- rma.mv(yi = lnRR_Ea, V = VCV_E2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_E3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -7.8702 15.7405 27.7405 42.6723 28.7649
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0109 0.1046 30 no Study_ID
## sigma^2.2 0.0319 0.1787 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 764.5984, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.4138, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_reinforcementAppetitive 0.2175 0.0726 2.9942 89 0.0036 0.0732
## Type_reinforcementAversive 0.2191 0.0410 5.3456 89 <.0001 0.1377
## Type_reinforcementNot applicable 0.0812 0.0560 1.4504 89 0.1505 -0.0300
## ci.ub
## Type_reinforcementAppetitive 0.3618 **
## Type_reinforcementAversive 0.3006 ***
## Type_reinforcementNot applicable 0.1924
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E3) ## R2_marginal R2_coditional
## 0.07720103 1.00000000
Reinforcement_E <-orchard_plot(mod_E3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34", "grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Reinforcement_E Fig. 4c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparisons
contra_mod_E3 <- contrast_fun(data = dat2, response = lnRR_Ea, moderator = Type_reinforcement,
VCV = VCV_E2)
res_table_mod_E3 <- get_estimate(model = mod_E3, contra = contra_mod_E3, moderator = Type_reinforcement)
res_table_mod_E3## # A tibble: 6 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Appetitive 0.217 0.0732 0.362 3.56e-3 -0.218 0.653
## 2 Aversive 0.219 0.138 0.301 6.89e-7 -0.200 0.638
## 3 Not applicable 0.0812 -0.0300 0.192 1.50e-1 -0.345 0.507
## 4 Appetitive-Aversive 0.00164 -0.163 0.167 9.84e-1 -0.442 0.445
## 5 Appetitive-Not applic… -0.136 -0.315 0.0425 1.33e-1 -0.585 0.312
## 6 Aversive-Not applicab… -0.138 -0.259 -0.0172 2.56e-2 -0.567 0.291
Age of enrichment
The age when individuals were exposed to environmental enrichment
dat3 <- filter(dat, Age_EE_exposure %in% c("Adult", "Adolescent"))
VCV_E3 <- impute_covariance_matrix(vi = dat3$lnRRV_E, cluster = dat3$Study_ID, r = 0.5)
mod_E4 <- rma.mv(yi = lnRR_Ea, V = VCV_E3, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat3)
summary(mod_E4) ##
## Multivariate Meta-Analysis Model (k = 88; method: REML)
##
## logLik Deviance AIC BIC AICc
## -6.9490 13.8980 23.8980 36.1698 24.6480
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0078 0.0885 29 no Study_ID
## sigma^2.2 0.0327 0.1808 88 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 86) = 782.6092, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 18.9060, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Age_EE_exposureAdolescent 0.1799 0.0370 4.8553 86 <.0001 0.1062 0.2535
## Age_EE_exposureAdult 0.2340 0.0620 3.7734 86 0.0003 0.1107 0.3573
##
## Age_EE_exposureAdolescent ***
## Age_EE_exposureAdult ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E4) ## R2_marginal R2_coditional
## 0.01127347 1.00000000
Age_E<- orchard_plot(mod_E4, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Age_E Fig. 4d Orchard plot showing the group-wise means of the categorical variable ‘Age_EE_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparisons
contra_mod_E4 <- contrast_fun(data = dat3, response = lnRR_Ea, moderator = Age_EE_exposure,
VCV = VCV_E3)
res_table_mod_E4 <- get_estimate(model = mod_E4, contra = contra_mod_E4, moderator = Age_EE_exposure)
res_table_mod_E4## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Adolescent 0.180 0.106 0.254 0.00000533 -0.227 0.587
## 2 Adult 0.234 0.111 0.357 0.000295 -0.185 0.653
## 3 Adolescent-Adult 0.0541 -0.0895 0.198 0.456 -0.371 0.479
Exercise enrichment
Did enrichment involve the addition of apparatus for voluntary exercise (e.g., a wheel or treadmill)?
mod_E5<- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~EE_exercise-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_E5)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -10.0303 20.0606 30.0606 42.5597 30.7749
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0096 0.0980 30 no Study_ID
## sigma^2.2 0.0345 0.1857 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 807.7169, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 16.1666, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_exerciseExercise 0.1849 0.0407 4.5464 90 <.0001 0.1041 0.2657
## EE_exerciseNo exercise 0.1900 0.0556 3.4151 90 0.0010 0.0795 0.3005
##
## EE_exerciseExercise ***
## EE_exerciseNo exercise ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E5) ## R2_marginal R2_coditional
## 0.0001360993 0.9999999987
Exercise_E <-orchard_plot(mod_E5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-0.5, 2) +
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Exercise_E Fig. 4d Orchard plot showing the group-wise means of the categorical variable ‘EE_exercise’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Pairwise comparisons
contra_mod_E5 <- contrast_fun(data = dat, response = lnRR_Ea, moderator = EE_exercise,
VCV = VCV_E)
res_table_mod_E5 <- get_estimate(model = mod_E5, contra = contra_mod_E5, moderator = EE_exercise)
res_table_mod_E5## # A tibble: 3 x 7
## Levels Estimate Lower_CI Upper_CI P_value Lower_PI Upper_PI
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Exercise 0.185 0.104 0.266 0.0000169 -0.240 0.610
## 2 No exercise 0.190 0.0795 0.301 0.000958 -0.241 0.621
## 3 Exercise-No exercise 0.00511 -0.132 0.142 0.941 -0.434 0.444
Publication bias & sensitivity analysis
Multi-moderator model
# filter data so that all K < 5 are removed
dat_Efm <- dat %>%
filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"), Type_reinforcement %in%
c("Appetitive", "Aversive", "Not applicable"), EE_social %in% c("Social",
"Non-social"), Age_EE_exposure %in% c("Adult", "Adolescent"))
VCV_Efm <- impute_covariance_matrix(vi = dat_Efm$lnRRV_E, cluster = dat_Efm$Study_ID,
r = 0.5)
mod_Efm <- rma.mv(yi = lnRR_Sa, V = VCV_Efm, mod = ~Type_assay - 1 + Learning_vs_memory +
Type_reinforcement + EE_social + EE_exercise + Age_EE_exposure, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat_Efm)
# summary(mod_Efm) r2_ml(mod_Efm)
res_Efm <- dredge(mod_Efm, trace = 2)
saveRDS(res_Efm, file = here("Rdata", "res_Efm.rds"))
# also saving the full model and data
saveRDS(mod_Efm, file = here("Rdata", "mod_Efm.rds"))
saveRDS(dat_Efm, file = here("Rdata", "dat_Efm.rds"))The akaike weights for the top set of models with AIC < 6
dat_Efm <- readRDS(file = here("Rdata", "dat_Efm.rds"))
mod_Efm <- readRDS(file = here("Rdata", "mod_Efm.rds"))
res_Efm <- readRDS(file = here("Rdata", "res_Efm.rds"))
res_Efm2 <- subset(res_Efm, delta <= 6, recalc.weights = FALSE)
importance(res_Efm2)## Age_EE_exposure Type_assay EE_exercise EE_social
## Sum of weights: 0.63 0.36 0.16 0.14
## N containing models: 11 10 5 5
## Learning_vs_memory Type_reinforcement
## Sum of weights: 0.10 0.07
## N containing models: 4 4
Funnel plot
Used to produce Fig. 7a
# funnel plot
Funnel_E <- funnel(mod_Efm, xlab = "lnRR", ylab = "Standard Error")# Funnel_EFig.7a Funnel plot showing the standard error and residuals (lnRR) from the full model.
# year published was scaled previously under stress PB
dat_Efm$sqrt_inv_e_n <- with(dat_Efm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
PB_MR_E <- rma.mv(lnRR_Sa, lnRRV_S, mods = ~1 + sqrt_inv_e_n + Learning_vs_memory +
Year_published + Type_assay + Type_reinforcement + EE_social + EE_exercise +
Age_stress_exposure, random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML",
test = "t", data = dat_Efm)
# PB_MR_E
estimates_PB_MR_E <- estimates.CI(PB_MR_E)
# estimates_PB_MR_ELeave-one-out analysis
Individually removing studies to distinguish that no singular study is having a disproportionate effect
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$Study_ID))) {
d <- dat %>%
filter(Study_ID != levels(dat$Study_ID)[i])
VCV_Eb <- impute_covariance_matrix(vi = d$lnRRV_E, cluster = d$Study_ID, r = 0.5)
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Ea, V = VCV_Eb, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", data = dat[dat$Study_ID !=
levels(dat$Study_ID)[i], ])
}
# TODO Shinichi to check
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_E0) {
df <- data.frame(est = mod_E0$b, lower = mod_E0$ci.lb, upper = mod_E0$ci.ub)
return(df)
}
# using dplyr to form data frame
MA_CVR_E <- lapply(LeaveOneOut_effectsize, function(x) est.func(x)) %>%
bind_rows %>%
mutate(left_out = levels(dat$Study_ID))
saveRDS(MA_CVR_E, file = here("Rdata", "MA_CVR_E.rds"))Used to produce Fig. S3
# telling ggplot to stop reordering factors
MA_CVR_E <- readRDS(file = here("Rdata", "MA_CVR_E.rds"))
MA_CVR_E$left_out <- as.factor(MA_CVR_E$left_out)
MA_CVR_E$left_out <- factor(MA_CVR_E$left_out, levels = MA_CVR_E$left_out)
# plotting
leaveoneout_E <- ggplot(MA_CVR_E) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_E0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_E0$b, lty = 1, lwd = 0.75, colour = "black") + geom_hline(yintercept = mod_E0$ci.ub,
lty = 3, lwd = 0.75, colour = "black") + geom_pointrange(aes(x = left_out, y = est,
ymin = lower, ymax = upper)) + xlab("Study left out") + ylab("lnRR, 95% CI") +
coord_flip() + theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
leaveoneout_Edat$Study_ID <- as.integer(dat$Study_ID)Fig. S3 Leave-one-group-out analysis showing meta-analytic mean and 95% CI when each individual study is removed from the data set.
Using SMD
Re-rerunning meta-analytic model using standardised mean difference (SMD) instead of lnRR
mod_E0a <- rma.mv(yi = SMD_Ea, V = VCV_Ea, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), test = "t", data = dat)
summary(mod_E0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -111.6284 223.2568 231.2568 241.3003 231.7220
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.2512 0.5012 30 no Study_ID
## sigma^2.2 0.4247 0.6517 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 673.4722, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.7241 0.1295 5.5895 91 <.0001 0.4668 0.9814 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E0a)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 8.673320e-01 3.223202e-01 5.450118e-01 1.749563e-09
Main effect: stress
Meta-analysis
mod_S0 <- rma.mv(yi = lnRR_Sa, V = VCV_S, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.1156 28.2312 36.2312 46.2747 36.6964
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0099 0.0993 30 no Study_ID
## sigma^2.2 0.0377 0.1941 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 946.9234, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.1052 0.0337 -3.1217 91 0.0024 -0.1721 -0.0383 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 9.376124e-01 1.946168e-01 7.429955e-01 1.975455e-10
orchard_plot(mod_S0, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) Fig. 3 Orchard plot showing meta-analytic mean and 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Meta-regression: uni-moderator
Learning vs Memory
Was the response learning or memory?
mod_S1 <- rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S1) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -14.5281 29.0562 39.0562 51.5553 39.7705
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0094 0.0970 30 no Study_ID
## sigma^2.2 0.0388 0.1969 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 946.3930, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 5.0145, p-val = 0.0086
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning -0.1211 0.0476 -2.5423 90 0.0127 -0.2157
## Learning_vs_memoryMemory -0.0974 0.0380 -2.5648 90 0.0120 -0.1728
## ci.ub
## Learning_vs_memoryLearning -0.0265 *
## Learning_vs_memoryMemory -0.0219 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S1) ## R2_marginal R2_coditional
## 0.002776034 1.000000000
LvsM_S <- orchard_plot(mod_S1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
LvsM_S Fig. 5a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Type of assay
The broad category of the type of assay used to measure learning or memory
dat$Type_assay<-as.factor(dat$Type_assay)
VCV_S1 <- impute_covariance_matrix(vi = dat1$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
mod_S2 <- rma.mv(yi = lnRR_Sa, V = VCV_S1, mod = ~Type_assay-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_S2)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -9.8028 19.6055 31.6055 46.5374 32.6299
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0161 0.1271 30 no Study_ID
## sigma^2.2 0.0279 0.1671 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 723.4973, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 6.7053, p-val = 0.0004
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Type_assayConditioning -0.0981 0.0375 -2.6192 89 0.0104 -0.1725 -0.0237
## Type_assayHabituation -0.4615 0.1126 -4.0969 89 <.0001 -0.6853 -0.2377
## Type_assayRecognition -0.0534 0.0645 -0.8287 89 0.4095 -0.1816 0.0747
##
## Type_assayConditioning *
## Type_assayHabituation ***
## Type_assayRecognition
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S2) ## R2_marginal R2_coditional
## 0.1853359 1.0000000
Learning_S <-orchard_plot(mod_S2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Learning_S Fig. 5b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Type of reinforcement
If conditioning was used, was aversive or appetitive reinforcement used?
VCV_S2 <- impute_covariance_matrix(vi = dat2$lnRRV_S, cluster = dat$Study_ID, r = 0.5)
mod_S3 <- rma.mv(yi = lnRR_Sa, V = VCV_S2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_S3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -13.8810 27.7621 39.7621 54.6939 40.7864
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0103 0.1016 30 no Study_ID
## sigma^2.2 0.0387 0.1966 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 920.8439, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.7942, p-val = 0.0130
##
## Model Results:
##
## estimate se tval df pval
## Type_reinforcementAppetitive -0.1846 0.0749 -2.4649 89 0.0156
## Type_reinforcementAversive -0.0730 0.0427 -1.7081 89 0.0911
## Type_reinforcementNot applicable -0.1172 0.0590 -1.9851 89 0.0502
## ci.lb ci.ub
## Type_reinforcementAppetitive -0.3334 -0.0358 *
## Type_reinforcementAversive -0.1579 0.0119 .
## Type_reinforcementNot applicable -0.2345 0.0001 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S3) ## R2_marginal R2_coditional
## 0.0366428 1.0000000
Reinforcement_S <-orchard_plot(mod_S3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Reinforcement_S Fig. 5c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Age of stress
The age when individuals were exposed to stress
mod_S4 <-rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Age_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_S4) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -12.4083 24.8166 38.8166 56.1579 40.2166
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0048 0.0694 30 no Study_ID
## sigma^2.2 0.0392 0.1979 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 881.1229, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 4.3703, p-val = 0.0029
##
## Model Results:
##
## estimate se tval df pval
## Age_stress_exposureAdolescent 0.0074 0.1159 0.0641 88 0.9490
## Age_stress_exposureAdult -0.2279 0.0622 -3.6664 88 0.0004
## Age_stress_exposureEarly postnatal -0.0561 0.0435 -1.2891 88 0.2007
## Age_stress_exposurePrenatal -0.1145 0.0743 -1.5404 88 0.1271
## ci.lb ci.ub
## Age_stress_exposureAdolescent -0.2228 0.2377
## Age_stress_exposureAdult -0.3514 -0.1044 ***
## Age_stress_exposureEarly postnatal -0.1425 0.0304
## Age_stress_exposurePrenatal -0.2621 0.0332
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S4) ## R2_marginal R2_coditional
## 0.09987307 1.00000000
Age_S <- orchard_plot(mod_S4, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Age_S Fig. 5d Orchard plot showing the group-wise means of the categorical variable ‘Age_stress_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Type of stress
The type of stressor used
dat5 <- filter(dat, Type_stress_exposure %in% c("Restraint", "Noise", "MS", "Combination"))
VCV_S3 <- impute_covariance_matrix(vi = dat5$lnRRV_S, cluster = dat5$Study_ID, r = 0.5)
mod_S4 <- rma.mv(yi = lnRR_Sa, V = VCV_S3, mod = ~Type_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat5)
summary(mod_S4) ##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -11.4138 22.8276 36.8276 53.5887 38.3618
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0115 0.1071 25 no Study_ID
## sigma^2.2 0.0401 0.2002 85 no ES_ID
## sigma^2.3 0.0000 0.0000 4 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 81) = 897.4553, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 81) = 2.8767, p-val = 0.0278
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_stress_exposureCombination -0.0500 0.0892 -0.5605 81 0.5767 -0.2274
## Type_stress_exposureMS -0.0539 0.0560 -0.9630 81 0.3384 -0.1653
## Type_stress_exposureNoise -0.1203 0.1036 -1.1608 81 0.2491 -0.3265
## Type_stress_exposureRestraint -0.2101 0.0704 -2.9863 81 0.0037 -0.3501
## ci.ub
## Type_stress_exposureCombination 0.1274
## Type_stress_exposureMS 0.0575
## Type_stress_exposureNoise 0.0859
## Type_stress_exposureRestraint -0.0701 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S4)## R2_marginal R2_coditional
## 0.07029554 1.00000000
Stressor<- orchard_plot(mod_S4, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Stressor Fig. 5e Orchard plot showing the group-wise means of the categorical variable ‘Type_stress_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Stess duration
Was the stress acute or chronic?
dat6 <- filter(dat, Stress_duration %in% c("Chronic", "Acute"))
VCV_S4 <- impute_covariance_matrix(vi = dat6$lnRRV_S, cluster = dat6$Study_ID, r = 0.5)
mod_S5 <-rma.mv(yi = lnRR_Sa, V = VCV_S4, mod = ~Stress_duration-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat6)
summary(mod_S5) ##
## Multivariate Meta-Analysis Model (k = 89; method: REML)
##
## logLik Deviance AIC BIC AICc
## -13.7898 27.5796 37.5796 49.9092 38.3204
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0062 0.0786 29 no Study_ID
## sigma^2.2 0.0391 0.1978 89 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 87) = 915.9393, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 87) = 6.6661, p-val = 0.0020
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Stress_durationAcute 0.0135 0.0659 0.2045 87 0.8384 -0.1176 0.1445
## Stress_durationChronic -0.1360 0.0373 -3.6456 87 0.0005 -0.2101 -0.0618
##
## Stress_durationAcute
## Stress_durationChronic ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S5) ## R2_marginal R2_coditional
## 0.08245896 1.00000000
Duration_S <- orchard_plot(mod_S5, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10))
Duration_S Fig. 5f Orchard plot showing the group-wise means of the categorical variable ‘Stress_duration’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Publication bias & sensitivity analysis
Multi-moderator model
# selecting moderator levels with k >=5
dat_Sfm <- dat %>%
filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"), Type_reinforcement %in%
c("Appetitive", "Aversive", "Not applicable"), Type_stress_exposure %in%
c("Restraint", "Noise", "MS", "Combination"), Stress_duration %in% c("Chronic",
"Acute"))
VCV_Sfm <- impute_covariance_matrix(vi = dat_Sfm$lnRRV_E, cluster = dat_Sfm$Study_ID,
r = 0.5)
mod_Sfm <- rma.mv(yi = lnRR_Sa, V = VCV_Sfm, mod = ~Type_assay - 1 + Learning_vs_memory +
Type_reinforcement + Type_stress_exposure + Age_stress_exposure + Stress_duration,
random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), test = "t", data = dat_Sfm)
# summary(mod_Sfm) r2_ml(mod_Sfm)
res_Sfm <- dredge(mod_Sfm, trace = 2)
saveRDS(res_Sfm, file = here("Rdata", "res_Sfm.rds"))
# also saving the full model and data
saveRDS(mod_Sfm, file = here("Rdata", "mod_Sfm.rds"))
saveRDS(dat_Sfm, file = here("Rdata", "dat_Sfm.rds"))The akaike weights for the top set of models with AIC < 6
dat_Sfm <- readRDS(file = here("Rdata", "dat_Sfm.rds"))
mod_Sfm <- readRDS(file = here("Rdata", "mod_Sfm.rds"))
res_Sfm <- readRDS(file = here("Rdata", "res_Sfm.rds"))
res_Sfm2 <- subset(res_Sfm, delta <= 6, recalc.weights = FALSE)
importance(res_Sfm2)## Type_assay Stress_duration Type_stress_exposure
## Sum of weights: 0.91 0.88 0.22
## N containing models: 7 6 2
## Learning_vs_memory Age_stress_exposure Type_reinforcement
## Sum of weights: 0.16 0.04 0.04
## N containing models: 2 1 1
Funnel plot
Used to produce Fig. 7b
# funnel plot
Funnel_S <- funnel(mod_Sfm, xlab = "lnRR", ylab = "Standard Error")Funnel_S## x y slab
## 1 -0.113053232 0.36468860 1
## 2 -0.015439898 0.20610340 2
## 3 -0.248010857 0.31117978 3
## 4 -0.232933972 0.39558646 4
## 5 -0.042387513 0.28415283 5
## 6 -0.371377055 0.24616892 6
## 7 -0.299243891 0.33372726 7
## 8 0.044834345 0.28792375 8
## 9 0.126615612 0.17056933 9
## 10 -0.177388895 0.16829598 10
## 11 -0.247450735 0.25159343 11
## 12 0.460979011 0.26162092 12
## 13 -0.554560458 0.33539070 13
## 14 0.680693862 0.30867691 14
## 15 -0.471382295 0.25172620 15
## 16 -0.394493983 0.23638600 16
## 17 -0.650018042 0.27849269 17
## 18 -0.300841465 0.25362984 18
## 19 0.027931562 0.19059106 19
## 20 0.080720040 0.19033477 20
## 21 0.041790076 0.24429542 21
## 22 0.203537585 0.10208024 22
## 23 0.124260874 0.14463442 23
## 24 0.116691372 0.14706024 24
## 25 -0.121449889 0.14833999 25
## 26 -0.001971635 0.19205933 26
## 27 -0.003907599 0.20022315 27
## 28 0.120523677 0.21113795 28
## 29 0.404659861 0.21388659 29
## 30 -0.109516461 0.22302661 30
## 31 -0.001061616 0.19288079 31
## 32 -0.004468007 0.19262936 32
## 33 0.043909789 0.18566567 33
## 34 0.115869661 0.18965184 34
## 35 -0.011278194 0.09467548 35
## 36 0.175061863 0.18048966 36
## 37 -0.122130273 0.22610402 37
## 38 0.002684314 0.17636045 38
## 39 0.035121843 0.18537686 39
## 40 -0.156416652 0.15860530 40
## 41 -0.138825592 0.15899389 41
## 42 0.053065175 0.17943678 42
## 43 0.138424978 0.18297543 43
## 44 0.102483572 0.18367269 44
## 45 0.170732200 0.19143563 45
## 46 0.139838423 0.18382247 46
## 47 0.128377093 0.17993531 47
## 48 0.221166287 0.18370137 48
## 49 0.044994292 0.18507026 49
## 50 -0.255771493 0.15305831 50
## 51 0.040084756 0.13800224 51
## 52 0.117345431 0.17133452 52
## 53 0.225493174 0.17017071 53
## 54 -0.041531740 0.18021370 54
## 55 0.380896989 0.79456429 55
## 56 0.014647194 0.18626025 56
## 57 -0.239220751 0.20755664 57
## 58 0.200890800 0.26208602 58
## 59 0.478075273 0.39791794 59
## 60 -0.693031272 0.44427913 60
## 61 -0.446007512 0.22443320 61
## 62 -0.385573815 0.31627073 62
## 63 0.187438535 0.29281093 63
## 64 -0.357529861 0.30767463 64
## 65 -0.206432699 0.18443973 65
## 66 -0.072514197 0.18259794 66
## 67 -0.058864366 0.39854817 67
## 68 0.552319003 0.28778615 68
## 69 0.499178177 0.26916891 69
## 70 0.527506647 0.32337519 70
## 71 -0.162795600 0.23927168 71
## 72 -0.032389680 0.19746325 72
## 73 0.107665252 0.24604039 73
## 74 0.209204201 0.16352434 74
## 75 -0.174302918 0.19626869 75
## 76 -0.021661369 0.19202699 76
## 77 -0.024341578 0.23180570 77
## 78 0.052275983 0.18465565 78
## 79 0.331124135 0.46142545 79
## 80 -0.102398694 0.20899450 80
## 81 0.349947062 0.56232826 81
## 82 0.203444436 0.16285272 82
Fig.7b Funnel plot showing the standard error and residuals (lnRR) from the full model.
# calculating inv effective sample size for use in full meta-regression
dat_Sfm$sqrt_inv_e_n <- with(dat_Sfm, sqrt(1/CC_n + 1/EC_n + 1/ES_n + 1/CS_n))
# time lag bias and eggers regression
dat_Sfm$c_Year_published <- as.vector(scale(dat_Sfm$Year_published, scale = F))
PB_MR_S <- rma.mv(lnRR_Sa, lnRRV_S, mods = ~1 + sqrt_inv_e_n + c_Year_published +
Type_assay + Learning_vs_memory + Type_reinforcement + Type_stress_exposure +
Age_stress_exposure, random = list(~1 | Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML",
test = "t", data = dat_Sfm, control = list(optimizer = "optim", optmethod = "Nelder-Mead"))
# PB_MR_S
estimates_PB_MR_S <- estimates.CI(PB_MR_S)
# estimates_PB_MR_SLeave-one-out sensitivity analysis
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$Study_ID))) {
d <- dat %>%
filter(Study_ID != levels(dat$Study_ID)[i])
VCV_Sb <- impute_covariance_matrix(vi = d$lnRRV_E, cluster = d$Study_ID, r = 0.5)
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_Sa, V = VCV_Sb, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", data = dat[dat$Study_ID !=
levels(dat$Study_ID)[i], ])
}
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_S0) {
df <- data.frame(est = mod_S0$b, lower = mod_S0$ci.lb, upper = mod_S0$ci.ub)
return(df)
}
# using dplyr to form data frame
MA_CVR_S <- lapply(LeaveOneOut_effectsize, function(x) est.func(x)) %>%
bind_rows %>%
mutate(left_out = levels(dat$Study_ID))
saveRDS(MA_CVR_S, file = here("Rdata", "MA_CVR_S.rds"))Used to produce Fig. S4
MA_CVR_S <- readRDS(file = here("Rdata", "MA_CVR_S.rds"))
# telling ggplot to stop reordering factors
MA_CVR_S$left_out <- as.factor(MA_CVR_S$left_out)
MA_CVR_S$left_out <- factor(MA_CVR_S$left_out, levels = MA_CVR_S$left_out)
# plotting
leaveoneout_S <- ggplot(MA_CVR_S) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_S0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_S0$b, lty = 1, lwd = 0.75, colour = "black") + geom_hline(yintercept = mod_S0$ci.ub,
lty = 3, lwd = 0.75, colour = "black") + geom_pointrange(aes(x = left_out, y = est,
ymin = lower, ymax = upper)) + xlab("Study left out") + ylab("lnRR, 95% CI") +
coord_flip() + theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
leaveoneout_Sdat$Study_ID <- as.integer(dat$Study_ID)Fig. S4 Leave-one-group-out analysis showing meta-analytic mean and 95% CI when each individual study is removed from the data set.
Using SMD
Re-rerunning meta-analytic model using standardised mean difference (SMD) instead of lnRR
mod_S0a <- rma.mv(yi = SMD_Sa, V = VCV_Sa, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), test = "t", data = dat)
summary(mod_S0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -120.2817 240.5634 248.5634 258.6068 249.0285
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.3874 0.6224 30 no Study_ID
## sigma^2.2 0.4969 0.7049 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 818.5592, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.4256 0.1498 -2.8403 91 0.0056 -0.7232 -0.1279 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S0a)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 8.959540e-01 3.924790e-01 5.034750e-01 1.438454e-09
Interaction: enrichment x stress
Meta-analysis
mod_ES0 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES0) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.8178 81.6355 89.6355 99.6790 90.1006
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0316 0.1777 30 no Study_ID
## sigma^2.2 0.0229 0.1513 92 no ES_ID
## sigma^2.3 0.0030 0.0544 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 303.2179, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1229 0.0596 2.0605 91 0.0422 0.0044 0.2414 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES0) ## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 0.81703063 0.44913873 0.32576306 0.04212884
orchard_plot(mod_ES0, mod = "Int", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 5, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 24), # change font sizes
legend.title = element_text(size = 15),
legend.text = element_text(size = 13)) Fig. 3 Orchard plot showing meta-analytic mean and 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Meta-regression: uni-moderator
Learning vs Memory
Was the response learning or memory?
mod_ES1 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES1) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.4769 80.9539 90.9539 103.4529 91.6682
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0292 0.1708 30 no Study_ID
## sigma^2.2 0.0232 0.1524 92 no ES_ID
## sigma^2.3 0.0034 0.0582 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 299.1854, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 2.9219, p-val = 0.0590
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Learning_vs_memoryLearning 0.1744 0.0722 2.4166 90 0.0177 0.0310
## Learning_vs_memoryMemory 0.1057 0.0619 1.7065 90 0.0914 -0.0174
## ci.ub
## Learning_vs_memoryLearning 0.3178 *
## Learning_vs_memoryMemory 0.2288 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES1) ## R2_marginal R2_coditional
## 0.0197648 0.9405080
LvsM_ES <- orchard_plot(mod_ES1, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
LvsM_ES Fig. 6a Orchard plot showing the group-wise means of the categorical variable ‘Learning_vs_memory’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Type of assay
The broad category of the type of assay used to measure learning or memory
VCV_ES1 <- impute_covariance_matrix(vi = dat1$lnRRV_ES, cluster = dat$Study_ID, r = 0.5)
mod_ES2 <- rma.mv(yi = lnRR_ESa, V = VCV_ES1, mod = ~Type_assay-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat1)
summary(mod_ES2)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0874 78.1748 90.1748 105.1066 91.1992
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0370 0.1924 30 no Study_ID
## sigma^2.2 0.0192 0.1386 92 no ES_ID
## sigma^2.3 0.0018 0.0422 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 293.9385, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.1062, p-val = 0.0305
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Type_assayConditioning 0.1525 0.0589 2.5877 89 0.0113 0.0354 0.2696
## Type_assayHabituation 0.1990 0.1415 1.4070 89 0.1629 -0.0820 0.4801
## Type_assayRecognition -0.0048 0.0800 -0.0606 89 0.9518 -0.1637 0.1541
##
## Type_assayConditioning *
## Type_assayHabituation
## Type_assayRecognition
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES2) ## R2_marginal R2_coditional
## 0.05775809 0.97105550
Learning_ES <- orchard_plot(mod_ES2, mod = "Type_assay", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Learning_ES Fig. 6b Orchard plot showing the group-wise means of the categorical variable ‘Type_assay’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Type of reinforcement
If conditioning was used, was aversive or appetitive reinforcement used?
VCV_ES2 <- impute_covariance_matrix(vi = dat2$lnRRV_ES, cluster = dat2$Study_ID, r = 0.5)
mod_ES3 <- rma.mv(yi = lnRR_ESa, V = VCV_ES2, mod = ~ Type_reinforcement-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat2)
summary(mod_ES3)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0604 78.1208 90.1208 105.0526 91.1452
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0382 0.1954 30 no Study_ID
## sigma^2.2 0.0189 0.1377 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 293.4724, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.2547, p-val = 0.0254
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_reinforcementAppetitive 0.1007 0.1075 0.9366 89 0.3515 -0.1129
## Type_reinforcementAversive 0.1573 0.0569 2.7618 89 0.0070 0.0441
## Type_reinforcementNot applicable 0.0147 0.0702 0.2101 89 0.8341 -0.1247
## ci.ub
## Type_reinforcementAppetitive 0.3143
## Type_reinforcementAversive 0.2704 **
## Type_reinforcementNot applicable 0.1541
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES3) ## R2_marginal R2_coditional
## 0.0586952 1.0000000
Reinforcement_ES <- orchard_plot(mod_ES3, mod = "Type_reinforcement", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Reinforcement_ES Fig. 6c Orchard plot showing the group-wise means of the categorical variable ‘Type_reinforcement’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Age of enrichment
The age when individuals were exposed to environmental enrichment
VCV_ES3 <- impute_covariance_matrix(vi = dat3$lnRRV_ES, cluster = dat3$Study_ID, r = 0.5)
mod_ES4 <- rma.mv(yi = lnRR_ESa, V = VCV_ES3, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat3)
summary(mod_ES4) ##
## Multivariate Meta-Analysis Model (k = 88; method: REML)
##
## logLik Deviance AIC BIC AICc
## -36.6407 73.2813 83.2813 95.5531 84.0313
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0332 0.1822 29 no Study_ID
## sigma^2.2 0.0207 0.1439 88 no ES_ID
## sigma^2.3 0.0007 0.0265 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 86) = 288.6493, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 86) = 3.1308, p-val = 0.0487
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Age_EE_exposureAdolescent 0.1284 0.0583 2.2006 86 0.0304 0.0124
## Age_EE_exposureAdult 0.1247 0.0921 1.3538 86 0.1793 -0.0584
## ci.ub
## Age_EE_exposureAdolescent 0.2444 *
## Age_EE_exposureAdult 0.3077
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES4) ## R2_marginal R2_coditional
## 0.0000400928 0.9871095822
Age_enrichment_ES <- orchard_plot(mod_ES4, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Age_enrichment_ES Fig. 6d Orchard plot showing the group-wise means of the categorical variable ‘Age_EE_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Exercise enrichment
Did enrichment involve the addition of apparatus for voluntary exercise (e.g., a wheel or treadmill)?
mod_ES5<- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~EE_exercise-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES5)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -40.4968 80.9935 90.9935 103.4926 91.7078
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0325 0.1803 30 no Study_ID
## sigma^2.2 0.0230 0.1517 92 no ES_ID
## sigma^2.3 0.0065 0.0805 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 90) = 297.3270, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 1.8401, p-val = 0.1647
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## EE_exerciseExercise 0.1051 0.0780 1.3474 90 0.1812 -0.0499 0.2602
## EE_exerciseNo exercise 0.1687 0.0967 1.7448 90 0.0844 -0.0234 0.3609
##
## EE_exerciseExercise
## EE_exerciseNo exercise .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES5) ## R2_marginal R2_coditional
## 0.01474459 0.89712469
Exercise_ES <- orchard_plot(mod_ES5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Exercise_ES Fig. 6d Orchard plot showing the group-wise means of the categorical variable ‘EE_exercise’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Age of stress
The age when individuals were exposed to stress
mod_ES7 <-rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Age_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES7) ##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0917 78.1834 92.1834 109.5247 93.5834
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0295 0.1718 30 no Study_ID
## sigma^2.2 0.0232 0.1522 92 no ES_ID
## sigma^2.3 0.0091 0.0954 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 88) = 286.4225, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 1.5932, p-val = 0.1832
##
## Model Results:
##
## estimate se tval df pval
## Age_stress_exposureAdolescent -0.0137 0.1762 -0.0775 88 0.9384
## Age_stress_exposureAdult 0.1677 0.1140 1.4708 88 0.1449
## Age_stress_exposureEarly postnatal 0.1067 0.0920 1.1602 88 0.2491
## Age_stress_exposurePrenatal 0.3179 0.1311 2.4254 88 0.0173
## ci.lb ci.ub
## Age_stress_exposureAdolescent -0.3639 0.3366
## Age_stress_exposureAdult -0.0589 0.3942
## Age_stress_exposureEarly postnatal -0.0761 0.2896
## Age_stress_exposurePrenatal 0.0574 0.5784 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES7) ## R2_marginal R2_coditional
## 0.09276574 0.86630830
Age_stress_ES<-orchard_plot(mod_ES7, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Age_stress_ES Fig. 6d Orchard plot showing the group-wise means of the categorical variable ‘Age_stress_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Type of stress
The type of stressor used
VCV_ES5 <- impute_covariance_matrix(vi = dat5$lnRRV_ES, cluster = dat5$Study_ID, r = 0.5)
mod_ES8 <- rma.mv(yi = lnRR_ESa, V = VCV_ES5, mod = ~Type_stress_exposure-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat5)
summary(mod_ES8)##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## logLik Deviance AIC BIC AICc
## -34.4046 68.8091 82.8091 99.5703 84.3434
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0426 0.2064 25 no Study_ID
## sigma^2.2 0.0232 0.1524 85 no ES_ID
## sigma^2.3 0.0104 0.1021 4 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 81) = 281.9708, p-val < .0001
##
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 81) = 0.5137, p-val = 0.7258
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Type_stress_exposureCombination 0.1111 0.1458 0.7618 81 0.4484 -0.1790
## Type_stress_exposureMS 0.1185 0.1099 1.0784 81 0.2840 -0.1001
## Type_stress_exposureNoise 0.1651 0.1795 0.9198 81 0.3604 -0.1920
## Type_stress_exposureRestraint 0.1374 0.1252 1.0978 81 0.2755 -0.1116
## ci.ub
## Type_stress_exposureCombination 0.4011
## Type_stress_exposureMS 0.3370
## Type_stress_exposureNoise 0.5221
## Type_stress_exposureRestraint 0.3865
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES8)## R2_marginal R2_coditional
## 0.004455703 0.863910284
Stressor_ES <- orchard_plot(mod_ES8, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34", "grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Stressor_ES Fig. 6e Orchard plot showing the group-wise means of the categorical variable ‘Type_stress_exposure’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Stress duration
Was the stress acute or chronic?
VCV_ES6 <- impute_covariance_matrix(vi = dat6$lnRRV_ES, cluster = dat6$Study_ID, r = 0.5)
mod_ES9 <-rma.mv(yi = lnRR_ESa, V = VCV_ES6, mod = ~Stress_duration-1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat6)
summary(mod_ES9) ##
## Multivariate Meta-Analysis Model (k = 89; method: REML)
##
## logLik Deviance AIC BIC AICc
## -35.9010 71.8020 81.8020 94.1315 82.5427
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0133 0.1155 29 no Study_ID
## sigma^2.2 0.0260 0.1611 89 no ES_ID
## sigma^2.3 0.0080 0.0895 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 87) = 278.4877, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 87) = 4.3877, p-val = 0.0153
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## Stress_durationAcute -0.0367 0.0930 -0.3946 87 0.6941 -0.2216 0.1482
## Stress_durationChronic 0.1854 0.0732 2.5347 87 0.0130 0.0400 0.3308
##
## Stress_durationAcute
## Stress_durationChronic *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES9) ## R2_marginal R2_coditional
## 0.1598311 0.8575918
Duration_ES<- orchard_plot(mod_ES9, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_colour_manual(values = c("grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34")) +
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Duration_ES Fig. 6f Orchard plot showing the group-wise means of the categorical variable ‘Stress_duration’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Order to treatment exposure
The order in which individuals were exposed to enrichment and stress
mod_ES10 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Exposure_order -1, random = list(~1|Study_ID,
~1|ES_ID,
~1|Strain),
test = "t",
data = dat)
summary(mod_ES10)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -39.0408 78.0817 90.0817 105.0135 91.1061
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0316 0.1777 30 no Study_ID
## sigma^2.2 0.0227 0.1507 92 no ES_ID
## sigma^2.3 0.0000 0.0000 6 no Strain
##
## Test for Residual Heterogeneity:
## QE(df = 89) = 292.2561, p-val < .0001
##
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.1546, p-val = 0.0287
##
## Model Results:
##
## estimate se tval df pval ci.lb
## Exposure_orderConcurrently 0.1492 0.1208 1.2351 89 0.2201 -0.0909
## Exposure_orderEnrichment first -0.1782 0.1659 -1.0744 89 0.2856 -0.5079
## Exposure_orderStress first 0.1370 0.0526 2.6046 89 0.0108 0.0325
## ci.ub
## Exposure_orderConcurrently 0.3893
## Exposure_orderEnrichment first 0.1514
## Exposure_orderStress first 0.2414 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES10)## R2_marginal R2_coditional
## 0.1027207 1.0000000
Order_ES <- orchard_plot(mod_ES10, mod = "Exposure_order", xlab = "lnRR", alpha=0.4) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
scale_colour_manual(values = c("grey34","grey34","grey34")) +
scale_fill_manual(values=c("grey34","grey34","grey34")) +
geom_point(aes(fill = name), size = 3, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
axis.text.x = element_text(size = 10), # change font sizes
axis.text.y = element_text(size = 10),
legend.title = element_text(size = 7),
legend.text = element_text(size = 7))
Order_ES Fig. 6g Orchard plot showing the group-wise means of the categorical variable ‘Exposure_order’ with their 95% confidence interval (think black line) and 95% prediction interval (thin black line). Individuals points show observed effect sizes and sample sizes.
Publication bias & sensitivity analysis
Multi-moderator model
dat_ESfm <- dat %>%
filter(Type_assay %in% c("Recognition", "Habituation", "Conditioning"), Type_reinforcement %in%
c("Appetitive", "Aversive", "Not applicable"), EE_social %in% c("Social",
"Non-social"), Age_EE_exposure %in% c("Adult", "Adolescent"), Type_stress_exposure %in%
c("Restraint", "Noise", "MS", "Combination"), Stress_duration %in% c("Chronic",
"Acute"), Age_stress_exposure %in% c("Prenatal", "Early postnatal", "Adult"))
VCV_ESfm <- impute_covariance_matrix(vi = dat_ESfm$lnRRV_ES, cluster = dat_ESfm$Study_ID,
r = 0.5)
mod_ESfm <- rma.mv(yi = lnRR_Sa, V = VCV_ESfm, mod = ~Type_assay - 1 + Learning_vs_memory +
Type_reinforcement + EE_social + EE_exercise + Age_EE_exposure + Type_stress_exposure +
Age_stress_exposure + Stress_duration + Exposure_order, random = list(~1 | Study_ID,
~1 | ES_ID, ~1 | Strain), test = "t", data = dat_ESfm)
# summary(mod_ESfm) r2_ml(mod_ESfm)
res_ESfm <- dredge(mod_ESfm, trace = 2)
saveRDS(res_ESfm, file = here("Rdata", "res_ESfm.rds"))
# also saving the full model and data
saveRDS(mod_ESfm, file = here("Rdata", "mod_ESfm.rds"))
saveRDS(dat_ESfm, file = here("Rdata", "dat_ESfm.rds"))The akaike weights for the top set of models with AIC < 6
dat_ESfm <- readRDS(file = here("Rdata", "dat_ESfm.rds"))
mod_ESfm <- readRDS(file = here("Rdata", "mod_ESfm.rds"))
res_ESfm <- readRDS(file = here("Rdata", "res_ESfm.rds"))
res_ESfm2 <- subset(res_ESfm, delta <= 6, recalc.weights = FALSE)
importance(res_ESfm2)## Type_assay Stress_duration Age_EE_exposure EE_exercise
## Sum of weights: 0.65 0.63 0.26 0.11
## N containing models: 19 18 7 6
## Learning_vs_memory EE_social Age_stress_exposure
## Sum of weights: 0.08 0.08 0.05
## N containing models: 4 4 2
## Type_stress_exposure Exposure_order
## Sum of weights: 0.03 0.02
## N containing models: 2 1
Funnel plot
Used to produce Fig. 7c
Funnel_ES <- funnel(mod_ESfm, xlab = "lnRR", ylab = "Standard Error")Funnel_ES## x y slab
## 1 -0.160575495 0.55343043 1
## 2 -0.274967420 0.70627539 2
## 3 -0.089909776 0.46527295 3
## 4 -0.418899318 0.36571674 4
## 5 -0.402782263 0.53501323 5
## 6 -0.058704027 0.45431531 6
## 7 0.098710844 0.15778981 7
## 8 -0.041805142 0.11742065 8
## 9 -0.234568128 0.44509952 9
## 10 0.473861619 0.42494140 10
## 11 -0.541677850 0.59839672 11
## 12 0.693576469 0.55379623 12
## 13 -0.506062429 0.32349206 13
## 14 -0.434662931 0.29560170 14
## 15 -0.752838102 0.41427250 15
## 16 -0.403661524 0.38640076 16
## 17 -0.074888497 0.14501213 17
## 18 -0.016611205 0.15220818 18
## 19 0.208016615 0.31484827 19
## 20 0.066418083 0.06641794 20
## 21 -0.016988105 0.13630983 21
## 22 0.102490150 0.27788492 22
## 23 0.077633012 0.12209543 23
## 24 0.265433711 0.18060215 24
## 25 0.549569896 0.20386929 25
## 26 0.029904759 0.23531678 26
## 27 0.036477985 0.10420416 27
## 28 0.033071595 0.10236053 28
## 29 0.154387233 0.13226659 29
## 30 0.068347398 0.10358055 30
## 31 -0.066651554 0.06766794 31
## 32 0.180810374 0.14954583 32
## 33 -0.116381761 0.30939379 33
## 34 0.008432825 0.12562534 34
## 35 0.040870354 0.17488050 35
## 36 0.039701546 0.09315773 36
## 37 0.119572535 0.09839159 37
## 38 0.083631129 0.10349806 38
## 39 0.151879757 0.14903289 39
## 40 0.120985979 0.10446433 40
## 41 0.031045847 0.10043162 41
## 42 0.118346228 0.11005125 42
## 43 0.155471736 0.12842298 43
## 44 -0.100968105 0.08874519 44
## 45 0.044022534 0.07229828 45
## 46 0.101114679 0.16689703 46
## 47 0.203773608 0.15471161 47
## 48 0.088563198 0.17334032 48
## 49 -0.005169248 0.17866865 49
## 50 0.434942304 0.36859082 50
## 51 0.554127069 0.73098123 51
## 52 -0.616979476 0.68762297 52
## 53 -0.364466901 0.25174124 53
## 54 -0.304033204 0.52455475 54
## 55 0.230466950 0.44507651 55
## 56 -0.314501446 0.47347958 56
## 57 -0.044976662 0.10592411 57
## 58 0.088941840 0.09196930 58
## 59 -0.100897815 0.73013998 59
## 60 0.510285555 0.59479999 60
## 61 0.451655914 0.34843255 61
## 62 0.479984384 0.70816567 62
## 63 -0.204829048 0.32322526 63
## 64 -0.079911943 0.14926781 64
## 65 0.065631803 0.33253908 65
## 66 0.186788247 0.09926520 66
## 67 -0.216336366 0.15728452 67
## 68 -0.069183632 0.12020637 68
## 69 -0.074226124 0.44173325 69
## 70 0.038912354 0.12603637 70
## 71 0.441601579 1.08275832 71
## 72 -0.149920957 0.21396887 72
## 73 0.302424799 1.01869334 73
## 74 0.205736931 0.09751054 74
Fig.7c Funnel plot showing the standard error and residuals (lnRR) from the full model.
Leave-one-out analysis
dat$Study_ID <- as.factor(dat$Study_ID)
LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$Study_ID))) {
d <- dat %>%
filter(Study_ID != levels(dat$Study_ID)[i])
VCV_ESb <- impute_covariance_matrix(vi = d$lnRRV_ES, cluster = d$Study_ID, r = 0.5)
LeaveOneOut_effectsize[[i]] <- rma.mv(yi = lnRR_ESa, V = VCV_ESb, random = list(~1 |
Study_ID, ~1 | ES_ID, ~1 | Strain), method = "REML", data = dat[dat$Study_ID !=
levels(dat$Study_ID)[i], ])
}
# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(mod_ES0) {
df <- data.frame(est = mod_ES0$b, lower = mod_ES0$ci.lb, upper = mod_ES0$ci.ub)
return(df)
}
# using dplyr to form data frame
MA_CVR_ES <- lapply(LeaveOneOut_effectsize, function(x) est.func(x)) %>%
bind_rows %>%
mutate(left_out = levels(dat$Study_ID))
saveRDS(MA_CVR_ES, , file = here("Rdata", "MA_CVR_ES.rds"))Used to produce Fig. S5
MA_CVR_ES <- readRDS(here("Rdata", "MA_CVR_ES.rds"))
# telling ggplot to stop reordering factors
MA_CVR_ES$left_out <- as.factor(MA_CVR_ES$left_out)
MA_CVR_ES$left_out <- factor(MA_CVR_ES$left_out, levels = MA_CVR_ES$left_out)
# plotting
leaveoneout_ES <- ggplot(MA_CVR_ES) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +
geom_hline(yintercept = mod_ES0$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
geom_hline(yintercept = mod_ES0$b, lty = 1, lwd = 0.75, colour = "black") + geom_hline(yintercept = mod_ES0$ci.ub,
lty = 3, lwd = 0.75, colour = "black") + geom_pointrange(aes(x = left_out, y = est,
ymin = lower, ymax = upper)) + xlab("Study left out") + ylab("lnRR, 95% CI") +
coord_flip() + theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))
leaveoneout_ESdat$Study_ID <- as.integer(dat$Study_ID)Fig. S5 Leave-one-group-out analysis showing meta-analytic mean and 95% CI when each individual study is removed from the data set.
Using SMD
Re-rerunning meta-analytic model using standardised mean difference (SMD) instead of lnR
mod_ES0a <- rma.mv(yi = SMD_ESa, V = VCV_ESa, random = list(~1 | Study_ID, ~1 | ES_ID,
~1 | Strain), test = "t", data = dat)
summary(mod_ES0a)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -126.0571 252.1141 260.1141 270.1576 260.5793
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.4653 0.6821 30 no Study_ID
## sigma^2.2 0.3698 0.6081 92 no ES_ID
## sigma^2.3 0.0000 0.0002 6 no Strain
##
## Test for Heterogeneity:
## Q(df = 91) = 257.4673, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.6880 0.1763 3.9017 91 0.0002 0.3377 1.0382 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES0a)## I2_total I2_Study_ID I2_ES_ID I2_Strain
## 6.657130e-01 3.709364e-01 2.947766e-01 3.351907e-08
‘Pairwise’ effect sizes
Enrichment relative to control
VCV_E20 <- impute_covariance_matrix(vi = dat$lnRRV_E2, cluster = dat$Study_ID, r = 0.5)
mod_E20 <- rma.mv(yi = lnRR_E2a, V = VCV_E20, random = list(~1 | Study_ID, ~1 | Strain,
~1 | ES_ID), test = "t", data = dat, control = list(optimizer = "optim", optmethod = "Nelder-Mead"))
summary(mod_E20)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -7.3235 14.6470 22.6470 32.6904 23.1121
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0037 0.0611 30 no Study_ID
## sigma^2.2 0.0000 0.0000 6 no Strain
## sigma^2.3 0.0281 0.1675 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 475.8327, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.1066 0.0291 3.6655 91 0.0004 0.0489 0.1644 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E20)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 8.610789e-01 1.011724e-01 2.607945e-09 7.599065e-01
orchard_plot(mod_E20, mod = "Int", xlab = "lnRR")Stress relative to control
VCV_S20 <- impute_covariance_matrix(vi = dat$lnRRV_S2, cluster = dat$Study_ID, r = 0.5)
mod_S20 <- rma.mv(yi = lnRR_S2a, V = VCV_S20, random = list(~1 | Study_ID, ~1 | Strain,
~1 | ES_ID), test = "t", data = dat)
summary(mod_S20)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -52.3561 104.7122 112.7122 122.7557 113.1773
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0323 0.1797 30 no Study_ID
## sigma^2.2 0.0000 0.0000 6 no Strain
## sigma^2.3 0.0798 0.2824 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 1003.0694, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.1846 0.0522 -3.5360 91 0.0006 -0.2882 -0.0809 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S20)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 9.489604e-01 2.734220e-01 9.879730e-10 6.755384e-01
orchard_plot(mod_S20, mod = "Int", xlab = "lnRR")Enrichment + stress relative to control
VCV_ES20 <- impute_covariance_matrix(vi = dat$lnRRV_ES2, cluster = dat$Study_ID,
r = 0.5)
mod_ES20 <- rma.mv(yi = lnRR_ES2a, V = VCV_ES20, random = list(~1 | Study_ID, ~1 |
Strain, ~1 | ES_ID), test = "t", data = dat)
summary(mod_ES20)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -10.3625 20.7250 28.7250 38.7684 29.1901
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0039 0.0625 30 no Study_ID
## sigma^2.2 0.0014 0.0377 6 no Strain
## sigma^2.3 0.0227 0.1508 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 402.2656, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.0801 0.0389 2.0594 91 0.0423 0.0028 0.1573 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES20)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 0.81513970 0.11355557 0.04132363 0.66026050
orchard_plot(mod_ES20, mod = "Int", xlab = "lnRR")Enrichment + stress relative to stress
VCV_E30 <- impute_covariance_matrix(vi = dat$lnRRV_E3, cluster = dat$Study_ID, r = 0.5)
mod_E30 <- rma.mv(yi = lnRR_E3a, V = VCV_E30, random = list(~1 | Study_ID, ~1 | Strain,
~1 | ES_ID), test = "t", data = dat)
summary(mod_E30)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -46.3447 92.6895 100.6895 110.7329 101.1546
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0235 0.1532 30 no Study_ID
## sigma^2.2 0.0263 0.1623 6 no Strain
## sigma^2.3 0.0543 0.2331 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 790.0249, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## 0.2465 0.1011 2.4389 91 0.0167 0.0457 0.4472 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E30)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 0.9456920 0.2132063 0.2391809 0.4933048
orchard_plot(mod_E30, mod = "Int", xlab = "lnRR")Enrichment + stress relative to enrichment
VCV_S30 <- impute_covariance_matrix(vi = dat$lnRRV_S3, cluster = dat$Study_ID, r = 0.5)
mod_S30 <- rma.mv(yi = lnRR_S3a, V = VCV_S30, random = list(~1 | Study_ID, ~1 | Strain,
~1 | ES_ID), test = "t", data = dat, control = list(optimizer = "optim", optmethod = "Nelder-Mead"))
summary(mod_S30)##
## Multivariate Meta-Analysis Model (k = 92; method: REML)
##
## logLik Deviance AIC BIC AICc
## -6.8788 13.7576 21.7576 31.8011 22.2228
##
## Variance Components:
##
## estim sqrt nlvls fixed factor
## sigma^2.1 0.0000 0.0000 30 no Study_ID
## sigma^2.2 0.0009 0.0292 6 no Strain
## sigma^2.3 0.0276 0.1662 92 no ES_ID
##
## Test for Heterogeneity:
## Q(df = 91) = 540.3522, p-val < .0001
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## -0.0087 0.0342 -0.2552 91 0.7992 -0.0766 0.0592
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S30)## I2_total I2_Study_ID I2_Strain I2_ES_ID
## 8.415428e-01 5.192043e-09 2.515933e-02 8.163835e-01
orchard_plot(mod_S30, mod = "Int", xlab = "lnRR")Figures
Panel of ‘focal’ ES and ‘pairwise’ ES orchard plots
Used to produce Fig. 3
mod_list1 <- list(mod_E0, mod_S0, mod_ES0)
mod_res1 <- lapply(mod_list1, function(x) mod_results(x, mod = "Int"))
merged1 <- submerge(mod_res1[[3]], mod_res1[[2]], mod_res1[[1]], mix = T)
merged1$mod_table$name <- factor(merged1$mod_table$name, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3"),
labels = rev(c("Enrichment ME", "Stress ME", "Interaction")))
merged1$data$moderator <- factor(merged1$data$moderator, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3"),
labels = rev(c("Enrichment ME", "Stress ME", "Interaction")))
orchard1<- orchard_plot(merged1, mod = "Int", xlab = "lnRR", angle = 0) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-2,4.5) +
geom_point(aes(fill = name), size = 4, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling +
scale_colour_manual(values = c("#00AEEF","#00A651","#ED1C24"))+ # change colours
scale_fill_manual(values=c("#00AEEF","#00A651","#ED1C24"))+
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10)) mod_list2 <- list(mod_S30, mod_E30, mod_ES20, mod_S20, mod_E20) #rearranged the order so that it matches intext results
mod_res2 <- lapply(mod_list2, function(x) mod_results(x, mod = "Int"))
merged2 <- submerge(mod_res2[[1]], mod_res2[[2]], mod_res2[[3]], mod_res2[[4]], mod_res2[[5]], mix = T)
merged2$mod_table$name <- factor(merged2$mod_table$name, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3", "Intrcpt4", "Intrcpt5"),
labels = rev(c("EC/CC", "CS/CC", "ES/CC", "ES/CS", "ES/EC")))
merged2$data$moderator <- factor(merged2$data$moderator, levels = c("Intrcpt1",
"Intrcpt2", "Intrcpt3", "Intrcpt4", "Intrcpt5"),
labels = rev(c("EC/CC", "CS/CC", "ES/CC", "ES/CS", "ES/EC")))
orchard2 <- orchard_plot(merged2, mod = "Int", xlab = "lnRR", angle = 0) +
geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
xlim(-2,4.5) +
geom_point(aes(fill = name), size = 4, shape = 21)+ # mean estimate
scale_size_continuous(range = c(1, 7))+ # change point scaling
scale_colour_manual(values = c("#7B81BE","#D7DF23","#F37158","#75CBF2","#97D2B4"))+ # change colours
scale_fill_manual(values=c("#7B81BE","#D7DF23","#F37158","#75CBF2","#97D2B4"))+
theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
text = element_text(size = 15), # change font sizes
legend.title = element_text(size = 10),
legend.text = element_text(size = 10)) p1 <- orchard1 + orchard2 + plot_annotation(tag_levels = "A")
p1# saved as PDF: 6 x 15 inchesFig. 3 Orchard plots showing the eight meta-analytic models. (A) main effects of enrichment and stress, and interaction of environmental enrichment and stress, (B) ‘pairwise’ comparisons between treatments and controls. Thick black lines = 95% CI, thin black lines = 95 % PI.
Panel of meta-regressions
Environmental enrichment
Used to produce Fig. 4
# Enrichment
E_mod <- (LvsM_E + Learning_E + Reinforcement_E)/(Age_E + Exercise_E + Social_E) +
plot_annotation(tag_levels = "A")
E_mod# saved as pdf 10 x 15 inchesFig. 4 Orchard plots showing the six different meta-regressions on the main effect of environmental enrichment on learning and memory. (A) learning versus memory response, (B) the type of assay, (C) the type of reinforcement, (D) the age at environmental enrichment, (E) type of manipulation of exercise during enrichment, (F) manipulation of the social environment during enrichment.
Stress
Used to produce Fig. 5
S_mod <- (LvsM_S + Learning_S + Reinforcement_S)/(Age_S + Stressor + Duration_S) +
plot_annotation(tag_levels = "A")
S_mod# saved as pdf 10 x 15 inchesFig. 5 Orchard plots showing the six different meta-regressions on the main effect of stress on learning and memory. (A) learning versus memory response, (B) the type of assay, (C) the type of reinforcement, (D) the age at stress, (E) the type of stressor (MS = maternal separation), (F) chronic or acute stress.
Interaction
Used to produce Fig. 6
ES_mod <- plot_grid(LvsM_ES, Learning_ES, Reinforcement_ES, Age_enrichment_ES, Age_stress_ES,
Order_ES, Exercise_ES, Social_ES, Stressor_ES, Duration_ES, labels = "AUTO",
ncol = 5)
ES_mod# saved as 10 x 20 inchesFig. 6 Orchard plots showing the 10 different meta-regressions of moderators on the interaction between environmental enrichment and stress learning and memory. (A) learning versus memory response, (B) the type of assay, (C) the type of reinforcement used, (D) the age at environmental enrichment, (E) the age at stress, (F) the order of treatment exposure, (G) if enrichment involved a manipulation of exercise, (H) manipulation of the social environment during enrichment, (I) the type of stressor, (F) stress was chronic or acute.
Panel of funnel plots
Used to produce Fig. 7
# EE
pdf(NULL)
dev.control(displaylist = "enable")
par(mar = c(4, 4, 0.1, 0))
A <- funnel(mod_Sfm, xlab = "Residuals (lnRR)", ylab = "Standard Error", xlim = c(-2,
2), ylim = c(0, 1.05))
A <- recordPlot()
invisible(dev.off())
# Stress
pdf(NULL)
dev.control(displaylist = "enable")
par(mar = c(4, 4, 0.1, 0))
B <- funnel(mod_Sfm, xlab = "Residuals (lnRR)", ylab = "Standard Error", xlim = c(-2,
2), ylim = c(0, 1.05))
B <- recordPlot()
invisible(dev.off())
# Interaction
pdf(NULL)
dev.control(displaylist = "enable")
par(mar = c(4, 4, 0.1, 0))
C <- funnel(mod_ESfm, xlab = "Residuals (lnRR)", ylab = "Standard Error", xlim = c(-2,
2), ylim = c(0, 1.05))
C <- recordPlot()
invisible(dev.off())
# putting together
ggdraw(A) + ggdraw(B) + ggdraw(C) + plot_annotation(tag_levels = "A")
# png(file = here('figs', 'Fig7_Funnels.png'))
# dev.off()knitr::include_graphics(here("figs", "funnels.png"))Fig. 7 Funnel plots of the standard error and residuals (lnRR) from the full models. (A) environmental enrichment main effect, (B) stress main effect, (C) environmental enrichment x stress interaction.
Software and package versions
sessionInfo() %>%
pander()R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
locale: en_AU.UTF-8||en_AU.UTF-8||en_AU.UTF-8||C||en_AU.UTF-8||en_AU.UTF-8
attached base packages: grid, stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: formatR(v.1.11), pander(v.0.6.4), gridGraphics(v.0.5-1), png(v.0.1-7), cowplot(v.1.1.1), ggthemr(v.1.1.0), ggalluvial(v.0.12.3), visdat(v.0.5.3), networkD3(v.0.4), GoodmanKruskal(v.0.0.3), patchwork(v.1.1.1), MuMIn(v.1.43.17), orchaRd(v.0.0.0.9000), clubSandwich(v.0.5.3), metafor(v.3.0-2), Matrix(v.1.3-3), here(v.1.0.1), forcats(v.0.5.1), stringr(v.1.4.0), dplyr(v.1.0.6), purrr(v.0.3.4), readr(v.1.4.0), tidyr(v.1.1.3), tibble(v.3.1.2), ggplot2(v.3.3.3) and tidyverse(v.1.3.1)
loaded via a namespace (and not attached): nlme(v.3.1-152), fs(v.1.5.0), lubridate(v.1.7.10), RColorBrewer(v.1.1-2), httr(v.1.4.2), rprojroot(v.2.0.2), tools(v.4.1.0), backports(v.1.2.1), utf8(v.1.2.1), R6(v.2.5.0), vipor(v.0.4.5), DBI(v.1.1.1), colorspace(v.2.0-1), withr(v.2.4.2), tidyselect(v.1.1.1), compiler(v.4.1.0), cli(v.2.5.0), rvest(v.1.0.0), pacman(v.0.5.1), xml2(v.1.3.2), sandwich(v.3.0-1), labeling(v.0.4.2), bookdown(v.0.22), scales(v.1.1.1), digest(v.0.6.27), rmarkdown(v.2.8), pkgconfig(v.2.0.3), htmltools(v.0.5.1.1), highr(v.0.9), dbplyr(v.2.1.1), htmlwidgets(v.1.5.3), rlang(v.0.4.11), readxl(v.1.3.1), rstudioapi(v.0.13), farver(v.2.1.0), generics(v.0.1.0), zoo(v.1.8-9), jsonlite(v.1.7.2), magrittr(v.2.0.1), ggbeeswarm(v.0.6.0), Rcpp(v.1.0.6), munsell(v.0.5.0), fansi(v.0.4.2), lifecycle(v.1.0.0), stringi(v.1.6.2), yaml(v.2.2.1), mathjaxr(v.1.4-0), crayon(v.1.4.1), lattice(v.0.20-44), haven(v.2.4.1), hms(v.1.1.0), knitr(v.1.33), pillar(v.1.6.1), igraph(v.1.2.6), codetools(v.0.2-18), stats4(v.4.1.0), reprex(v.2.0.0), glue(v.1.4.2), evaluate(v.0.14), modelr(v.0.1.8), vctrs(v.0.3.8), rmdformats(v.1.0.2), cellranger(v.1.1.0), gtable(v.0.3.0), assertthat(v.0.2.1), xfun(v.0.23), broom(v.0.7.6), beeswarm(v.0.4.0) and ellipsis(v.0.3.2)
Social enrichment
Did enrichment involve more conspecifics (note that we did not include studies that only provided social enrichment but no other form of abiotic enrichment)?
Pairwise comparisons